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ABSTRACT 

A predictor-corrector algorithm applicable to the numerical 
integration of an artificial earth satellite in multirevolution steps 
is reviewed. The multirevolution predictor and corrector formulas 
are developed to yield general recursive relations for the coeffi- 
cients, which are similar to those for the Adams’ integration 
formulas. In addition, the applicability of the algorithm to orbit 
calculations is demonstrated numerically. 
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NUMERICAL INTEGRATION OF ORBITS IN 
MULTIREVOLUTION STEPS 


INTRODUCTION 

Many artificial satellite orbit determination problems, such as lifetime 
studies or computing long-range predictions, require the numerical integration 
of the equations of motion over many periods of revolution. An algorithm de- 
signed to save computing time in solving such problems has been proposed in a 
paper by Cohen and Hubbard (1960). Essentially, the method involves combining 
the numerical integration with a procedure which "steps" the calculations ahead 
in multirevolution increments. This "stepping" procedure, both in formulation 
and usage, is similar to the Adams’ predictor-corrector process in numerical 
integration. Briefly, a cycle in the algorithm consists of first extrapolating or 
"predicting" the orbital elements n-revolutions ahead. Then, starting with these 
extrapolated values, the equations of motion are integrated over one revolution. 
Finally, the extrapolated values are improved by using a corrector formula along 
with the information obtained from the single revolution integration. 

The multirevolution "predictor" or extrapolation formula was introduced in 
a paper by Mace and Thomas (1960). In this paper it was demonstrated that if 
successive revolutions are similar, so that the successive descending nodes are 
close together, one may extrapolate without loss of accuracy as many revolutions 
as will make the change in the orbital elements from node to extrapolated node 
about the same as the change in step by step integration over one revolution. 
Supplementing this work, Cohen and Hubbard developed a corrector formula for 
this n-revolution extrapolation. In addition, it was also shown that this corrector 
formula reduces to the Adams' formula in the limiting case, i.e., the coefficients 
of the corrector were shown to be polynomials in 1/n with the constant terms 
being the Adams’ coefficients. 

In this paper, the similarity in derivation between the multirevolution for- 
mulas and the Adams' predictor-corrector formulas will be explored and used 
to develop recursive formulas for the coefficients. In addition, the multirevo- 
lution algorithm will be detailed and examples demonstrating its usage will be 
given. 


FORMULATION 

Let f be a function defined on an interval [x Q , x m J with equal spacing h (i.e., 
x. - x. . = h, i = 1, 2, • • • m) and assume f is continuously differentiable. 


1 


Consider the difference operators V, A, E a , D, 1 defined by 


Vf(x) - f(x) - f(x-h) (backward difference) 


Af(x) - f(x+h) - f(x) (forward difference) 


E a f(x) = f ( x + ah ) 


(shifting operator) 


Df(x) = f' (x) 


(differentiating operator) 


1 f(x) = f(x) 


( identity) 


where xe [x 0 , x m J . Also, if L is any one of these operators, define L n by 

L n f(x) = L(L n-1 f(x)). 


Some basic relations between these operators are* 


A = 

V 

(1) 

1 - v 

E = 

1 + A 

(2) 

E = 

e hD . 

(3) 


We begin by deriving the Adams’ integration formulas. The predictor can 
be derived by expressing the forward difference operator in terms of backward 


•Equality interpreted in the sense of operators. For verification of these formulas and justification for the 
algebraic manipulation of these operators see Hildebrand, pp. 128-134. 
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differences of the derivative. We have from (2) and (3) that 


hD - - log ( 1 - V) 


(4) 


and hence, using (1), 


V V 1 

A - 1 - V - 1 - V - log(l -V) hD 




D 


^i=0 


(5) 


where the coefficients a. can be determined as follows: 


Using the expressions 


1 - x 


x 


00 

Z> 


i=0 


(6) 


and 


00 

\ 1 x 1 

- log ( 1 “ x ) = X / r^TT 


i=0 


(7) 


we have 


V 1 

1 - V - log ( 1 - V) 


00 

L 


i=0 


V* 


' 1 v 1 


i=0 


00 

= z>‘ vl 


i=0 
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where the a. are given recursively by 


a o ~ 1 


a . 


- £> ‘-i 
j=l 


> (8) 


i = 1, 2, 3, 


where 

= ro • 

Hence, using (5) and (8), our formula emerges in the form 


( 9 ) 


Af (x) = f ( x + h) - f(x) 


= h 




V+ T2 y2 + •••} f ' ( x ) - 


( 10 ) 


which is the familiar Adams -Bashforth formula used to obtain a "predicted" 
value of f (x + h) when solving 1 st order initial value problems. The corre- 
sponding corrector formula can be derived by expressing the backward differ- 
ence operator in terms of backward differences of the derivative: 

As above, using (4), we have 


V 


V 

- log ( 1 - V) 



where the coefficients a.* are given recursively by 


a* 


1 


i 




(ii) 


(12) 


where the /3. are given by (9). 
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We have then, using (11) and (12), 

Vf(x) - f(x) - f(x-h) 

= h {»- 1 V - T2 V2 - •••}'■ <*> • ( 13 ) 

which is the Adams-Moulton formula used to obtain a corrected value of f(x) in 
the numerical integration process. 

We note that from (5) and (11) we have 


[l + V + V 2 + • • • J £a 0 * + a* V + • • • J - a Q + a.j V + • • • 
so equating coefficients of like powers, we have the relationship 



(14) 


We now wish to derive the multirevolution predictor and corrector formulas. 
This is done precisely as above, except that differences of the derivative are re- 
placed by differences at a large interval of differences at a small interval. To 
this end, we define the differences A and V by 

n n 


A n f ( x ) ' f ( x + nh) " f(x) 


V n f O) = f ( x ) " f(x -nh) , 


where n is same fixed integer, i.e., A n and V n are the forward and backward 
differences taken at n times the step of A and V. A relation between V n and A is 
then given by 


(l-V„) = (1 +A)“ n 
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since 


(l- V „)f(x) = f(x -nh) 


= E" n f ( x ) 

= (1 + A)" n f(x) . 


So we have, in analogy with (4), 

A = (! -V n)' 1/n ■ 1 - (15) 

and derive the multirevolution predictor or extrapolation formula by expressing 
the forward difference A in terms of backward differences A 1 of the forward 

n n 

difference A. We have from (1) and (15), in analogy with (5), 



(16) 


where the coefficients y. , which are polynomials in l/n, can be determined using 
series (6) and the expansion 


( 1 - x) a = (i) xi • 

i=0 


We have 


(i-v„)- ,/ " - 1 = 2Z(-1V CY n )v„‘ 


(17) 





i=0 
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where b = 1 and 


and hence 


IK-h-0 


C" 1 )" ( k + 1)! 


k = 1,2, 3, 




n (1-Vj-Vn-i 


00 

n V Vi 
/ , n 
i = 0 

co 


where the y. are given recursively by 


= 1 


■ £ b i 


Using (16) and (19), we then obtain 

A n f ( x ) = f ( x +nh ) " f (x) 


00 

n J2 7i ^ 


i = 1,2, 3, 


{i + 4(i-K)v„ + n(s-|^)v„ 2 


+ ' “r nA f(x) , (20) 
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which is the formula used to compute the "predicted" or extrapolated value of 
f ( x +nh) by using the backward differences V n * of the forward difference of f(x) 
at l/n th the step. 

Again, in analogy with the integration formulas, the corresponding multi- 
revolution corrector can be derived by expressing the backward difference V n in 
terms of backward differences V 1 of A. 

n 

We have (see (11)), 


V 


V 


(l -V n ) 1/n -l 


A 




/►A 


( 21 ) 


.i=0 


where, as above, we find that the coefficients y* are given recursively by 


y o* = 1 




i 

?i* = " 21 ^ 7 *i 

i = i 


( 22 ) 


i = 1,2, 3, 


where the b. are given by (18). 
Using (21) and (22) we obtain 


V n f O) - f ( x ) " f(x -nh) 





V 2 _ . . . 


•jnAf(x) , ( 23 ) 


which is the formula used to compute a corrected value of f(x) by using a for- 
ward difference of an extrapolated value for f(x). 
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We note the similarity between the coefficients given by (8), (19) and (12), 
(22). We also note, in analogy with (13), the relation 


k 

7 > = 7 k • 

i=0 


Finally, we remark that in practice, an algebraic equivalent of Equations (10) 
and (13) known as the "summed" form of the difference formulae is often used. 

It has been established (Henrici, 1962) that this modification considerably re- 
duces the propagated round-off error associated with these formulae. Formally, 
one can obtain this modification by applying the inverse operator V -1 to both 
sides of Equations (10) and (13). In like manner, one can obtain the "summed" 
form of the extrapolation formulae by applying the operator V n “ 1 to both sides of 
Equations (20) and (23). 


THE MULTIREVOLUTION ALGORITHM 

Let f . denote the value of an orbital element at the descending node of the 
j th revolution and let n be the number of revolutions to be "stepped." Also let 
k be the order of the highest backward difference V 1 to be retained in formulas 
(16) and (21). We then have the multirevolution formulas: 


A 

n 


f . 

i 



v 

n 


f . 
j 



(24) 


(25) 


We see that, in analogy with the Adams' integration process, a set of start- 
ing values must be computed before these formulas can be used. Assume then 
that the values Af = f , + 1 - f . , i = 0, n, 2n, 3n, • • • , kn, have been computed, 
i.e., we assume knowledge of the orbital elements at descending nodes 0, 1, n, 
n + 1, 2n, 2n + 1, • • • kn, kn + 1. (See table on next page). Using these values, 
we can compute the differences 

V „‘ K „) • i = 1. 2, ••• k . 
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v „ ) 

= Af. 

kn 

- Af 

(k-l)n 

VK.) 

= 

kn 

- 2 Af + Af 

* *(k-l)n 1 (k-2)n 


etc. We can then form a "starting" table (entries above the diagonal): 



Figure 1 
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The multirevolution cycle is then: 

(a) Extrapolate the orbital elements n revolutions ahead using (24) with 

j = nk, i.e., we compute a "predicted” value for f the element at the de- 

scending node of the (k + l)n th revolution. 

(b) Starting with these extrapolated values, we integrate one revolution to 
obtain the elements at descending node (k + l)n + 1 and form the forward difference 


Af 


(k + l)n 


f 


(k +l)n+l 


f 


(k +l)n ’ 


(c) We form the backward differences (Af ) , i = 1, 2, • • • k and use 
formula (25) with j = (k + 1 )n to obtain a corrected value for f , and the 
cycle is complete. The process is then repeated with j = (k + l)n in (a). 

We make the following remarks concerning the computation procedure: 

(i) Values of the orbital elements at the descending node of any revolution 
can be obtained by inverse interpolation, i.e., first, a direct interpolation is 
made to find the time of nodal crossing, followed by an inverse interpolation to 
obtain the values of the orbital elements. 

(ii) Other points in the orbit such as perigee, apogee or ascending node 
could be used instead of the descending node. 

(iii) The time associated with the extrapolated orbital elements can also be 
obtained by the extrapolation process, allowing one to obtain values for the orbital 
elements at any point in a computed revolution. 

(iv) For any particular orbit, the number of revolutions which can be 
"stepped” by this process and yet maintain some predictable accuracy can be 
determined by comparing the change in the orbital elements from node to node 
with the change required by the short step integration. 

(v) For any particular orbit, the number of differences V^ 1 ) to be retained 
in the extrapolation formulas can be determined precisely as in the case of the 
integration formula, i.e., retaining differences of sufficiently high order so that 
the "truncated" terms are within the accuracy of the computations. 


NUMERICAL RESULTS 

The equations of motion considered in the numerical computations are of the 
form 

^ ux — — 

X - - + Fj + F 2 , (26) 
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where x = (x, y, z), R = (x 2 +y 2 + z 2 ) 1/2 , Fj and F 2 perturbation functions, and 
fj. is a constant. The formulation of these perturbations and the values of the re- 
lated constants which were used are given in the appendix. 

Before beginning, we make the following remarks concerning the numerical 
results: 

(i) All short-step numerical integrations were performed with an Adams- 
Cowell 13 th order predictor-corrector process (Henrici, pp. 192-194, 291-293), 
using a Runge-Kutta type starter, so that this process was used both for the 
comparison integrations (Table A), and for the short-step integrations required 
by the multirevolution algorithm. 

(ii) A predictor-corrector tolerance of 1 x 10" 12 was used for all the 
Cowell integrations. The number of derivative evaluations tabulated in the tables 
below was obtained by multiplying the number of steps taken at the Cowell time 
step h by the average number of predictor-corrector iterations taken over the 
range of integration. 

(iii) For the examples, three orbits were used where in each case the in- 
tegrations were performed for approximately 100 revolutions: 


Initial Elements* 

Example I 

Example II 

Example III 

Semi-major axis (earth radii) 

1.15 

1.26 

6.71 

Eccentricity 

.075 

.072 

.003 

Inclination 

1.52 

1.03 

.0004 

Mean anomaly 

6.03 

3.71 

3.80 

Arg. of perigee 

1.15 

3.14 

0.31 

Long, of ascending node 

4.76 

6.16 

2.29 

Period (mins.) 

104.25 

120.37 

1469.46 

Range of integration (mins.) 

11,000 

12,500 

148,000 


(iv) All computations were performed on the UNIVAC 1108 computer using 
double-precision arithmetic. 

In Table A, the results of the comparison integrations are tabulated. Error 
estimates were obtained by integrating the Equations (26) with F t = F 2 = 0, and 
comparing with the Keplerian solution. For the cases considered, the approximate 
ratio of the change in the position which occurred at the Cowell time step h, and 
the corresponding change which occurred from descending node to descending 
node, was found to be 10 or greater. 


•Angles in radians. 
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Table A 


Orbit 

h 

(mins.) 

No. of Rev. 
Computed 

Total No. of 
Der. Eval. 

Estimated 

Error 

Computation Time 
(mins.) 

I 

1.0 

105 

17,337 

2 x 10 -9 

1.14 

II 

.8 

104 

15,617 

4 x 10‘ 12 

1.09 

III 

22.0 

100 

7,257 

9 x 10" 9 

.71 


We now consider the results obtained by applying the multirevolution algorithm 
to our test cases. The elements used in the extrapolation process were the po- 
sition and velocity vectors themselves. (The time associated with the extrapo- 
lated values was also obtained in this manner.) We begin by tabulating the results 
obtained by using the predictor formula (24) only (Table B). The error estimate 
given is the difference between the position at the node of the 100 th revolution 
obtained by the comparison integration (Table A) and by the multirevolution 
process. In each case, for each n, the order of the process (k) was increased 
until no change in accuracy occurred. 


Table B 


Orbit 

n 

k 

Total No. of 
Der. Eval. 

Estimated 

Error 

Computation Time 
(mins.) 

I 

3 

2 

6049 

1 x lO -8 

.55 



4 

6737 

3 x 10" 9 

.59 


5 

2 

4586 

3 x 10” 8 

.39 



4 

5934 

2 x 10" 9 

.48 



2 

4337 

8 x 10' 8 

.35 



4 

6343 

4 x 10~ 9 

OO 

• 



2 

4542 

1 x 10" 7 

.36 



4 

7209 

1 x 10 -9 

.52 

II 


2 

5650 

9 x 10" 8 

.51 



4 

6268 

8 x 10" 11 

.55 



2 

4263 

5 x 10' 7 

.36 



4 

5483 

6 x 10”“ 

.44 



2 

4013 

1 x 10' 6 

.33 



4 

5834 

4 x 10' 10 

.45 



6 

7656 

8 x 10~ n 

.57 


9 

2 

4188 

3 x 10' 6 

.33 



4 

6611 

2 x 10' 9 

.49 



6 

9034 

4 x 10' 10 

.65 

III 


2 

2329 

6 x 10' 7 

.62 



2 

1803 

5 x 10' 7 

.45 



2 

1815 

5 x 10' 7 

.40 


9 

2 

1923 

5 x 10' 7 

.37 
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Examining these results, we find that significant gains in efficiency were 
obtained by application of the algorithm to the numerical integrations, with no 
significant loss in accuracy. 

We remark that a good portion of the computation time tabulated was used in 
forming the starting table (Figure 1), and in restarting the short-step integrations 
required by the process, since these were done with the Runge-Kutta, Adams- 
Cowell process. 

In Table C, the results obtained by applying the predictor (24), followed by 
an application of the corrector (25), are tabulated. As in Table B, for each n, the 
order (k) was increased until no change in accuracy occurred. 


Table C 


Orbit 

n 

a 

Total No. of 
Der. Eval. 

Estimated 

Error 

Computation Time 
(mins.) 

I 

3 

2 

10,903 

2 x 10" 9 

.99 


5 

2 

7,317 

5 x 10" 9 

.65 



4 

8,361 

3 x 10' 9 

.70 


n 

2 

6,157 

1 x 10“ 8 

.52 


1 

4 

7,860 

1 x 10' 9 

.62 


I 

2 

5,907 

2 x 10' 8 

OO 

• 


1 

4 

8,271 

2 x 10" 9 

.62 

II 

1 

2 

10,194 

2 x 10' 8 

.92 


1 

4 

10,528 

8 x 10" 11 

.94 



2 

6,819 

7 x 10“ 8 

.59 


1 1 

4 

7,755 

7xl0' n 

.65 



2 

5,717 

2 x 10 -7 

OO 

• 



4 

7,254 

6 x 10' n 

.58 


9 

2 

5,466 

4 x 10 -7 

.44 



4 

7,605 

1 x 10 _1 ° 

.58 



6 

9,744 

2 x 10‘ n 

.72 

in 


2 

4,086 

6 x 10‘ 7 

.99 



2 

2,781 

5 x 10" 7 

.68 



2 

2,442 

5 x 10" 7 

.56 


D 

2 

2,387 

5 x 10" 7 

.48 


Examining these results, we find that an application of the corrector yielded 
slight gains in accuracy, but at some expense in efficiency. Again, as in Table B, 
this cost in efficiency was due to the process used for all the short-step inte- 
gration required by the algorithm. 
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At this point, several considerations should be made: 

(1) As suggested in References 1 and 4, an iterative process using the multi- 
revolution formulas could be used to form the required starting table (Figure 1), 
instead of the short-step integration process. Such a scheme would increase the 
efficiency of both the predictor and predictor-corrector processes used to obtain 
the results. 

(2) Since much of the computation time, especially in Table C, was used in 
restarting the short-step integrations by the Runge-Kutta technique an increase 
in efficiency could be obtained by extrapolating the required values for these 
integrations, as well as the nodal points. 

(3) With more efficient starting processes, it could be feasible to iterate on 
the corrector (25) to improve the accuracy when using large n. This would be 
analogous to the Adams' integration process. 


SUMMARY AND CONCLUSIONS 

The multirevolution algorithm has been reviewed and found to be an effective 
tool applicable to the orbit determination problem. 

The recursive formulas obtained for the coefficients required by the process 
have been found to be an effective and simple means of obtaining these coefficients 
for any given n and k. The numerical results have indicated that the "predictor 
only" scheme can be applied to the common multi-step integration process with 
a Runge-Kutta starter to yield considerable gains in efficiency and, if more ef- 
ficient starting processes were used, gains in accuracy due to an application of 
the corrector could be obtained at little, or no expense in efficiency. 
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APPENDIX 


THE EQUATIONS OF MOTION 


The equations of motion considered in this paper are of the form 


x ( t o) = x o 
^o) = *o 

where x = (x, y, z), R = | x| = ( x 2 + y 2 + z 2 ) 1/2 , i± is a constant, x Q , x Q the 
given initial position and velocity vectors, F x the perturbation due to the non- 
sphericity of the earth, and F 2 the perturbation due to drag. 


is given by 

F x = ^ [jR 2 (5 s-l) + Hz(7s-3) + | (-63s 2 +42s-3)] 

F y = ^ [JR 2 (5s- 1) + Hz(7s-3) + | (-63s 2 +42s-3)] 

F z = [j R 2 (5s - 1) + ^ (- 63s 2 + 70s- 15)] + (35s 2 - 30s + 3 ) 


t'Rt 


V ' if* 
i 1 \ O 


page blank not filmed. 
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where s - (z/R)J and J, H, K are the 2 nd , 3 rd and 4 th harmonics of the earth’s 
potential field. F 2 is of the form 

c d A 

2M H V r | V r 

where V r is the relative velocity of the satellite, p the atmospheric density* at 
the satellite position, and A, M, C D are the cross-sectional area, mass, and drag 
coefficient of the satellite respectively. 

Numerical constants: 

Unit of Time = 806.832 secs. 

Unit of Length = 6378.388 kms. 

M = 1 

J = +0.162 x 10" 2 
H = -0.640 x 10" 5 
K = +0.690 x 10" s 
A = +0.116 x 10 s cm 2 
M = +0.158 x 10 6 gm 
C D = +2.3 


♦The density was computed (from a standard model) for satellite heights of less than 1000 kms; otherwise p 
was set to zero. Hence, in our examples, only orbit I was influenced by F 2 * 
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